Using the package glmmTMB to model groups of species individually.
| month | site | salinity_ppt | salinity_regime | quadrat_no | balanus_pt | chthamalus_pt | barnacle_recruits_pt | barnacles_total_pt | mytilus_pt | fucus_pt | ulva_pt | hildenbrandia_pt | mastocarpus_pt | mastocarpus_crust_pt | pyropia_pt | endocladia_pt | microcladia_pt | polysiphonia_pt | diatoms_pt | urospora_pt | acrosiphonia_pt | melanosiphon_pt | leathesia_pt | greens_pt | greens_pyropia_pt | reds_pt | reds_fucus_pt | littorina_spp_no | sitkana_no | littorina_total_no | paradigitalis_no | unknown_limpet_no | pelta_no | scutum_no | persona_no | digitalis_no | all_limpets_no | gastropods_no | anemone_pt | hemigrapsis_no | hermit_crab_no |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| May | Ruckle | 30 | H | 1 | 1.00 | 0.50 | 0.0 | 1.50 | 0.250 | 0 | 0.00 | 0.00 | 0.00 | 0.50 | 0.00 | 0 | 0.00 | 0 | 0 | 0 | 0 | 0 | 0 | 0.00 | 0.00 | 0.50 | 0.50 | 41 | 0 | 41 | 49 | 0 | 0 | 0 | 0 | 2 | 51 | 92 | 0 | 0 | 0 |
| May | Ruckle | 30 | H | 2 | 8.25 | 13.00 | 0.0 | 21.25 | 0.025 | 0 | 0.00 | 0.50 | 0.00 | 0.00 | 0.00 | 0 | 0.00 | 0 | 0 | 0 | 0 | 0 | 0 | 0.00 | 0.00 | 0.50 | 0.50 | 98 | 0 | 98 | 22 | 0 | 1 | 0 | 0 | 0 | 23 | 121 | 0 | 0 | 0 |
| May | Ruckle | 30 | H | 3 | 1.50 | 26.50 | 0.0 | 28.00 | 0.250 | 0 | 0.00 | 1.00 | 0.00 | 2.50 | 0.75 | 0 | 0.50 | 0 | 0 | 0 | 0 | 0 | 0 | 0.00 | 0.75 | 4.75 | 4.75 | 66 | 0 | 66 | 32 | 0 | 1 | 0 | 0 | 0 | 33 | 99 | 0 | 0 | 0 |
| May | Ruckle | 30 | H | 4 | 4.50 | 18.75 | 18.5 | 41.75 | 0.025 | 0 | 0.00 | 0.00 | 0.25 | 3.00 | 7.25 | 0 | 0.00 | 0 | 0 | 0 | 0 | 0 | 0 | 0.00 | 7.25 | 10.50 | 10.50 | 58 | 0 | 58 | 13 | 0 | 0 | 0 | 0 | 1 | 14 | 72 | 0 | 0 | 0 |
| May | Ruckle | 30 | H | 5 | 1.50 | 1.25 | 0.0 | 2.75 | 0.500 | 0 | 0.00 | 0.25 | 0.00 | 0.75 | 71.50 | 0 | 0.00 | 0 | 0 | 0 | 0 | 0 | 0 | 0.00 | 71.50 | 72.50 | 72.50 | 0 | 0 | 0 | 15 | 0 | 0 | 0 | 0 | 1 | 16 | 16 | 5 | 0 | 0 |
| May | Ruckle | 30 | H | 6 | 8.25 | 11.00 | 0.0 | 19.25 | 0.025 | 0 | 0.25 | 0.00 | 0.00 | 0.25 | 15.00 | 0 | 0.25 | 0 | 0 | 0 | 0 | 0 | 0 | 0.25 | 15.25 | 15.50 | 15.50 | 6 | 0 | 6 | 16 | 0 | 1 | 0 | 0 | 0 | 17 | 23 | 9 | 0 | 0 |
Models fit to the data:
1. Zero-inflated Poisson model with a single zero-inflation parameter by applying to all observations (ziformula~1)
2. Zero-inflated Negative Binomial model with NB2 parameterization (variance = \(\mu\)(1 + \(\mu\)/\(\kappa\)), increasing quadratically with the mean)
3. Zero-inflated Negative Binomial model with NB1 parameterization (variance = \(\mu\) + \(\phi\)\(\mu\), increasing linearly with the mean)
4. Poisson model
5. Tweedie model
6. Gaussian model
7. Negative Binomial with NB2 parameterization
8. Negative Binomial with NB1 parameterization
Questions
Do I need to include salinity regime twice in the model syntax? i.e. salinity_regime + (1|salinity_regime/site) YES
Abundance of different limpet species in survey dataset
Fitting the models, and using AIC to select best model for the data:
limpet_zipoisson <- glmmTMB(all_limpets_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
limpet_ziNB <- update(limpet_zipoisson, family = nbinom2)
limpet_ziNB1 <- update(limpet_zipoisson, family = nbinom1) # best model
limpet_poisson <- glmmTMB(all_limpets_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
limpet_tweedie <- glmmTMB(all_limpets_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
limpet_gaussian <- glmmTMB((all_limpets_no + 1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
limpet_NB <- update(limpet_poisson, family = nbinom2)
limpet_NB1 <- update(limpet_poisson, family = nbinom1)
Model Comparison
| dAIC | df | |
|---|---|---|
| limpet_ziNB1 | 0.00000 | 12 |
| limpet_ziNB | 4.31957 | 12 |
| limpet_NB | 5.76179 | 11 |
| limpet_NB1 | 12.68349 | 11 |
| limpet_tweedie | 16.21041 | 12 |
| limpet_gaussian | 420.95556 | 11 |
| limpet_zipoisson | 1971.27929 | 11 |
| limpet_poisson | 2357.48378 | 10 |
Adding a dispersion formula for salinity regime to account for heterogeneity
limpet_ziNB1_2 <- glmmTMB(all_limpets_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, dispformula = ~salinity_regime, family = nbinom1)
#checking assumptions
limpet_sim_res_2 <- simulateResiduals(limpet_ziNB1_2)
plotResiduals(limpet_sim_res_2, form = survey_data$salinity_regime) # good
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 2.3297650 | 0.2711593 | 8.5918671 | 0.0000000 |
| monthJune | -0.1177717 | 0.2678733 | -0.4396546 | 0.6601873 |
| monthJuly | -0.5723872 | 0.3116878 | -1.8364117 | 0.0662968 |
| monthSeptember | -1.7281043 | 0.3864513 | -4.4717257 | 0.0000078 |
| salinity_regimeH | 1.2788265 | 0.3512972 | 3.6402978 | 0.0002723 |
| monthJune:salinity_regimeH | 0.2112509 | 0.3231194 | 0.6537858 | 0.5132498 |
| monthJuly:salinity_regimeH | 0.5471796 | 0.3623420 | 1.5101193 | 0.1310130 |
| monthSeptember:salinity_regimeH | 2.1470353 | 0.4238716 | 5.0652966 | 0.0000004 |
| term | statistic | df | p.value |
|---|---|---|---|
| (Intercept) | 73.82018 | 1 | 0.0000000 |
| month | 22.21336 | 3 | 0.0000589 |
| salinity_regime | 13.25177 | 1 | 0.0002723 |
| month:salinity_regime | 27.51028 | 3 | 0.0000046 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | 0.1177717 | 0.2678733 | 179 | 0.4396546 | 0.9715086 |
| May - July | L | 0.5723872 | 0.3116878 | 179 | 1.8364117 | 0.2597615 |
| May - September | L | 1.7281043 | 0.3864513 | 179 | 4.4717257 | 0.0000806 |
| June - July | L | 0.4546154 | 0.3156435 | 179 | 1.4402813 | 0.4760002 |
| June - September | L | 1.6103326 | 0.3907160 | 179 | 4.1214907 | 0.0003330 |
| July - September | L | 1.1557171 | 0.4194462 | 179 | 2.7553405 | 0.0325067 |
| May - June | H | -0.0934792 | 0.1807377 | 179 | -0.5172091 | 0.9548836 |
| May - July | H | 0.0252075 | 0.1853411 | 179 | 0.1360061 | 0.9990990 |
| May - September | H | -0.4189310 | 0.1741913 | 179 | -2.4050059 | 0.0797715 |
| June - July | H | 0.1186867 | 0.1823445 | 179 | 0.6508927 | 0.9151248 |
| June - September | H | -0.3254519 | 0.1706535 | 179 | -1.9070912 | 0.2287283 |
| July - September | H | -0.4441386 | 0.1747989 | 179 | -2.5408543 | 0.0571591 |
limpet abundance from transects in high and low salinity sites over the summer
littorina spp, sitkana, and littorina total
I’m using total littorines since the maximum number of sitkana is much lower and it’s mostly just Horseshoe Bay that has sitkana.
littorine_zipoisson <- glmmTMB(littorina_total_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
littorine_ziNB <- update(littorine_zipoisson, family = nbinom2)
littorine_ziNB1 <- update(littorine_zipoisson, family = nbinom1)
littorine_poisson <- glmmTMB(littorina_total_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
littorine_gaussian <- glmmTMB((littorina_total_no+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
littorine_tweedie <- update(littorine_poisson, family = tweedie)
littorine_NB <- update(littorine_poisson, family = nbinom2)
littorine_NB1 <- update(littorine_poisson, family = nbinom1)
Littorine model selection using AIC. The best model is the one using a negative binomial distribution, without modelling the zeroes. I needed to add a dispersion formula to model the variance in each salinity region separately.
## dAIC df
## littorine_NB1 0.0 11
## littorine_ziNB1 2.0 12
## littorine_NB 22.7 11
## littorine_ziNB 23.5 12
## littorine_tweedie 67.7 12
## littorine_gaussian 872.5 11
## littorine_zipoisson 11274.6 11
## littorine_poisson 14479.6 10
Adding a dispersion formula to account for heterogeneity in the high and low salinity regions.
# adding in a dispersion formula
littorine_NB1_2 <- glmmTMB(littorina_total_no ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, dispformula = ~salinity_regime*month, family = nbinom1)
#checking assumptions
set.seed(10)
littorine_sim_res_2 <- simulateResiduals(littorine_NB1_2)
par(mfrow = c(1,2))
plotResiduals(littorine_sim_res_2, form = survey_data$salinity_regime) # good
plotResiduals(littorine_sim_res_2, form = survey_data$month) # good
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 3.9557704 | 0.5896000 | 6.7092439 | 0.0000000 |
| monthJune | -0.2261032 | 0.6831038 | -0.3309939 | 0.7406491 |
| monthJuly | -2.8960237 | 0.6007762 | -4.8204703 | 0.0000014 |
| monthSeptember | -1.0691721 | 0.7068567 | -1.5125725 | 0.1303883 |
| salinity_regimeH | 0.0223592 | 0.7374873 | 0.0303181 | 0.9758133 |
| monthJune:salinity_regimeH | -0.5203907 | 0.7339034 | -0.7090725 | 0.4782795 |
| monthJuly:salinity_regimeH | 2.8659941 | 0.6523130 | 4.3935874 | 0.0000111 |
| monthSeptember:salinity_regimeH | 0.9981662 | 0.7608822 | 1.3118537 | 0.1895695 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 9.902932 | 3 | 0.0194095 |
| salinity_regime | 2.380059 | 1 | 0.1228929 |
| month:salinity_regime | 27.946482 | 3 | 0.0000037 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | 0.2261032 | 0.6831038 | 174 | 0.3309939 | 0.9874522 |
| May - July | L | 2.8960237 | 0.6007762 | 174 | 4.8204703 | 0.0000184 |
| May - September | L | 1.0691721 | 0.7068567 | 174 | 1.5125725 | 0.4321789 |
| June - July | L | 2.6699205 | 0.6935247 | 174 | 3.8497843 | 0.0009466 |
| June - September | L | 0.8430689 | 0.7872045 | 174 | 1.0709655 | 0.7076616 |
| July - September | L | -1.8268516 | 0.7169325 | 174 | -2.5481502 | 0.0562093 |
| May - June | H | 0.7464939 | 0.2682973 | 174 | 2.7823381 | 0.0302530 |
| May - July | H | 0.0300296 | 0.2541264 | 174 | 0.1181679 | 0.9994080 |
| May - September | H | 0.0710059 | 0.2815942 | 174 | 0.2521569 | 0.9943548 |
| June - July | H | -0.7164643 | 0.2498075 | 174 | -2.8680663 | 0.0237927 |
| June - September | H | -0.6754880 | 0.2777027 | 174 | -2.4324139 | 0.0748013 |
| July - September | H | 0.0409763 | 0.2640371 | 174 | 0.1551915 | 0.9986641 |
Littorine abundance from transect data in high and low salinity sites over the summer
species: Balanus, Chthalamus, recruits and total
For percent cover data, I’m rounding to nearest integer.
Best model is negative binomial distribution
balanus_poisson <- glmmTMB(round(balanus_pt) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
balanus_NB <- update(balanus_poisson, family = nbinom2)
balanus_NB1 <- update(balanus_poisson, family = nbinom1) #best model
balanus_gaussian <- glmmTMB((balanus_pt+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
balanus_tweedie <- glmmTMB(balanus_pt ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
# models with zero-inflation
balanus_zipoisson <- glmmTMB(round(balanus_pt) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
balanus_ziNB <- update(balanus_zipoisson, family = nbinom2)
balanus_ziNB1 <- update(balanus_zipoisson, family = nbinom1)
AICtab(balanus_poisson, balanus_NB, balanus_NB1, balanus_gaussian, balanus_tweedie, balanus_zipoisson, balanus_ziNB, balanus_ziNB1)
## dAIC df
## balanus_NB1 0.0 11
## balanus_ziNB1 2.0 12
## balanus_NB 35.4 11
## balanus_ziNB 37.4 12
## balanus_tweedie 44.5 12
## balanus_gaussian 170.4 11
## balanus_zipoisson 1554.2 11
## balanus_poisson 1817.0 10
Checking model assumptions:
#checking assumptions
set.seed(6789)
balanus_sim_res <- simulateResiduals(balanus_NB1)
plot(balanus_sim_res)
# plot residuals against time
par(mfrow = c(1,2))
plotResiduals(balanus_sim_res, form = survey_data$month) # different levels of variance in each month
#plot residuals against salinity regime
plotResiduals(balanus_sim_res, form = survey_data$salinity_regime) # more variance in the low salinity sites
Adding a dispersion formula for salinity regime and month
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.2848, p-value = 0.6995
## alternative hypothesis: true autocorrelation is not 0
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 1.9913985 | 0.3883563 | 5.1277619 | 0.0000003 |
| monthJune | 0.4464807 | 0.3207715 | 1.3918962 | 0.1639538 |
| monthJuly | 0.8889510 | 0.2854156 | 3.1145847 | 0.0018420 |
| monthSeptember | 0.7842295 | 0.3862107 | 2.0305743 | 0.0422982 |
| salinity_regimeH | 0.3468004 | 0.5304337 | 0.6538054 | 0.5132372 |
| monthJune:salinity_regimeH | -0.0164852 | 0.4184566 | -0.0393952 | 0.9685753 |
| monthJuly:salinity_regimeH | 0.0815283 | 0.3724967 | 0.2188700 | 0.8267513 |
| monthSeptember:salinity_regimeH | 0.2481843 | 0.4724092 | 0.5253588 | 0.5993338 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 32.4966405 | 3 | 0.0000004 |
| salinity_regime | 0.8598830 | 1 | 0.3537714 |
| month:salinity_regime | 0.4067149 | 3 | 0.9388519 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | -0.4464807 | 0.3207715 | 174 | -1.3918962 | 0.5061165 |
| May - July | L | -0.8889510 | 0.2854156 | 174 | -3.1145847 | 0.0114717 |
| May - September | L | -0.7842295 | 0.3862107 | 174 | -2.0305743 | 0.1807860 |
| June - July | L | -0.4424703 | 0.2654686 | 174 | -1.6667520 | 0.3444467 |
| June - September | L | -0.3377489 | 0.3717126 | 174 | -0.9086291 | 0.8002707 |
| July - September | L | 0.1047215 | 0.3416693 | 174 | 0.3064994 | 0.9899779 |
| May - June | H | -0.4299955 | 0.2687221 | 174 | -1.6001492 | 0.3812266 |
| May - July | H | -0.9704793 | 0.2393569 | 174 | -4.0545283 | 0.0004367 |
| May - September | H | -1.0324139 | 0.2720511 | 174 | -3.7949269 | 0.0011582 |
| June - July | H | -0.5404839 | 0.2128022 | 174 | -2.5398419 | 0.0574004 |
| June - September | H | -0.6024184 | 0.2490077 | 174 | -2.4192765 | 0.0772000 |
| July - September | H | -0.0619345 | 0.2169907 | 174 | -0.2854248 | 0.9918680 |
Balanus abundance (% cover) in from transects in high and low salinity sites over the summer
Only 1 plot at Rope Site has Chthamalus in it, and no plots at Lions Bay. This is causing convergence issues
chthamalus_poisson <- glmmTMB(round(chthamalus_pt) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log")) #convergence problems
chthamalus_NB <- update(chthamalus_poisson, family = nbinom2)
chthamalus_NB1 <- update(chthamalus_poisson, family = nbinom1) # convergence issues: This often occurs when a model is overparameterized (i.e. the data does not contain information to estimate the parameters).
chthamalus_gaussian <- glmmTMB((chthamalus_pt+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
chthamalus_tweedie <- glmmTMB(chthamalus_pt ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log")) #convergence issues
# models with zero-inflation
chthamalus_zipoisson <- glmmTMB(round(chthamalus_pt) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log")) #convergence issues
chthamalus_ziNB <- glmmTMB(round(chthamalus_pt) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = nbinom2(link = "log")) # best model
chthamalus_ziNB1 <- update(chthamalus_ziNB, family = nbinom1) #convergence issues
AICtab(chthamalus_NB, chthamalus_gaussian, chthamalus_ziNB)
## dAIC df
## chthamalus_ziNB 0.0 12
## chthamalus_NB 1.3 11
## chthamalus_gaussian 563.6 11
Checking model assumptions:
Adding a dispersion formula for salinity regime.
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.8696123 | 0.7872713 | -1.1045904 | 0.2693371 |
| monthJune | -19.5969277 | 4449.3512293 | -0.0044044 | 0.9964858 |
| monthJuly | -0.6649044 | 0.9308482 | -0.7142995 | 0.4750420 |
| monthSeptember | -0.4624549 | 0.9193908 | -0.5030014 | 0.6149633 |
| salinity_regimeH | 3.8193920 | 0.9228073 | 4.1388835 | 0.0000349 |
| monthJune:salinity_regimeH | 19.9082169 | 4449.3512355 | 0.0044744 | 0.9964300 |
| monthJuly:salinity_regimeH | 0.5627986 | 0.9603736 | 0.5860205 | 0.5578617 |
| monthSeptember:salinity_regimeH | 0.2304269 | 0.9509023 | 0.2423245 | 0.8085287 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 6.1112070 | 3 | 0.1063233 |
| salinity_regime | 22.7604991 | 1 | 0.0000018 |
| month:salinity_regime | 0.3606038 | 3 | 0.9482550 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | 19.5969277 | 4449.3512293 | 179 | 0.0044044 | 1.0000000 |
| May - July | L | 0.6649044 | 0.9308482 | 179 | 0.7142995 | 0.8913227 |
| May - September | L | 0.4624549 | 0.9193908 | 179 | 0.5030014 | 0.9582756 |
| June - July | L | -18.9320233 | 4449.3512080 | 179 | -0.0042550 | 1.0000000 |
| June - September | L | -19.1344729 | 4449.3512047 | 179 | -0.0043005 | 1.0000000 |
| July - September | L | -0.2024495 | 0.8285204 | 179 | -0.2443507 | 0.9948560 |
| May - June | H | -0.3112892 | 0.2332882 | 179 | -1.3343543 | 0.5423276 |
| May - July | H | 0.1021057 | 0.2352839 | 179 | 0.4339682 | 0.9725488 |
| May - September | H | 0.2320279 | 0.2418300 | 179 | 0.9594674 | 0.7725801 |
| June - July | H | 0.4133949 | 0.2302449 | 179 | 1.7954576 | 0.2789042 |
| June - September | H | 0.5433171 | 0.2350804 | 179 | 2.3111970 | 0.0992972 |
| July - September | H | 0.1299222 | 0.2369264 | 179 | 0.5483653 | 0.9468890 |
Chthamalus abundance from transect data in low and high salinity sites over the summer
Species:
Pyropia, Hildenbrandia, Mastocarpus, Mastocarpus crust, Endocladia, Microcladia, and Polysiphonia
Notes
Mostly Mastocarpus and Pyropia
Not very much red in the Sept/low plots.
Looks like there isn’t a lot at Rope Site.
red_poisson <- glmmTMB(round(reds) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
red_NB <- update(red_poisson, family = nbinom2)
red_NB1 <- update(red_poisson, family = nbinom1) #best model
red_gaussian <- glmmTMB((reds+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
red_tweedie <- glmmTMB(reds ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
# models with zero-inflation
red_zipoisson <- glmmTMB(round(reds) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
red_ziNB <- update(red_zipoisson, family = nbinom2) #convergence issues
red_ziNB1 <- update(red_zipoisson, family = nbinom1)
There is basically no red algae at Rope Site, which is making it very difficult to fit a model that meets the assumption of homogeneity of variance.
## dAIC df
## red_NB1 0.0 11
## red_ziNB1 2.0 12
## red_NB 19.6 11
## red_tweedie 78.9 12
## red_zipoisson 508.9 11
## red_gaussian 509.1 11
## red_poisson 792.1 10
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.1046, p-value = 0.2805
## alternative hypothesis: true autocorrelation is not 0
Adding a dispersion parameter for salinity_regime:
Adding dispersion parameter for site instead:
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 1.3955417 | 0.3310490 | 4.2155134 | 0.0000249 |
| monthJune | -0.5682088 | 0.4505585 | -1.2611212 | 0.2072652 |
| monthJuly | -0.3527778 | 0.4088760 | -0.8627989 | 0.3882481 |
| monthSeptember | -0.9547281 | 0.5077820 | -1.8801931 | 0.0600818 |
| salinity_regimeH | 0.6056305 | 0.3783527 | 1.6007037 | 0.1094426 |
| monthJune:salinity_regimeH | 0.7737543 | 0.4948812 | 1.5635154 | 0.1179314 |
| monthJuly:salinity_regimeH | 0.4175858 | 0.4610854 | 0.9056582 | 0.3651168 |
| monthSeptember:salinity_regimeH | 0.6503174 | 0.5589100 | 1.1635457 | 0.2446082 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 6.490347 | 3 | 0.0900440 |
| salinity_regime | 12.655822 | 1 | 0.0003744 |
| month:salinity_regime | 2.850151 | 3 | 0.4153106 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | 0.5682088 | 0.4505585 | 176 | 1.2611212 | 0.5888808 |
| May - July | L | 0.3527778 | 0.4088760 | 176 | 0.8627989 | 0.8239858 |
| May - September | L | 0.9547281 | 0.5077820 | 176 | 1.8801931 | 0.2402994 |
| June - July | L | -0.2154311 | 0.4594790 | 176 | -0.4688595 | 0.9657853 |
| June - September | L | 0.3865193 | 0.5523964 | 176 | 0.6997136 | 0.8970712 |
| July - September | L | 0.6019503 | 0.5177022 | 176 | 1.1627347 | 0.6511123 |
| May - June | H | -0.2055455 | 0.2047147 | 176 | -1.0040583 | 0.7472381 |
| May - July | H | -0.0648080 | 0.2131860 | 176 | -0.3039975 | 0.9902161 |
| May - September | H | 0.3044107 | 0.2335373 | 176 | 1.3034778 | 0.5619437 |
| June - July | H | 0.1407375 | 0.2003558 | 176 | 0.7024377 | 0.8960095 |
| June - September | H | 0.5099562 | 0.2225810 | 176 | 2.2911036 | 0.1040051 |
| July - September | H | 0.3692187 | 0.2288667 | 176 | 1.6132477 | 0.3738178 |
Red Algae (percent cover) from transect data in low and high salinity sites over the summer
Species: ulva_pt, urospora_pt, acrosiphonia_pt
Notes
Mostly Ulva
green_poisson <- glmmTMB(round(greens) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
green_NB <- update(green_poisson, family = nbinom2)
green_NB1 <- update(green_poisson, family = nbinom1) #best model
green_gaussian <- glmmTMB((greens+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
green_tweedie <- glmmTMB(greens ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
# models with zero-inflation
green_zipoisson <- glmmTMB(round(greens) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
green_ziNB <- update(green_zipoisson, family = nbinom2)
green_ziNB1 <- update(green_zipoisson, family = nbinom1) #convergence issues
Checking model assumptions:
## dAIC df
## green_NB1 0.0 11
## green_NB 11.3 11
## green_ziNB 13.3 12
## green_tweedie 88.0 12
## green_zipoisson 462.5 11
## green_gaussian 610.4 11
## green_poisson 951.5 10
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.1239, p-value = 0.8461
## alternative hypothesis: true autocorrelation is not 0
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 0.3910871 | 0.5332449 | 0.7334101 | 0.4633084 |
| monthJune | 0.2256104 | 0.4656709 | 0.4844845 | 0.6280420 |
| monthJuly | 0.8711570 | 0.4315772 | 2.0185428 | 0.0435348 |
| monthSeptember | 1.9120327 | 0.4043576 | 4.7285685 | 0.0000023 |
| salinity_regimeH | 0.3165282 | 0.7172861 | 0.4412858 | 0.6590061 |
| monthJune:salinity_regimeH | -0.2971707 | 0.6931032 | -0.4287539 | 0.6681023 |
| monthJuly:salinity_regimeH | -0.7566436 | 0.6465476 | -1.1702829 | 0.2418871 |
| monthSeptember:salinity_regimeH | -0.9512348 | 0.5879712 | -1.6178253 | 0.1057002 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 39.3806196 | 3 | 0.0000000 |
| salinity_regime | 0.3592205 | 1 | 0.5489395 |
| month:salinity_regime | 3.2021357 | 3 | 0.3614974 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | -0.2256104 | 0.4656709 | 181 | -0.4844845 | 0.9624619 |
| May - July | L | -0.8711570 | 0.4315772 | 181 | -2.0185428 | 0.1849684 |
| May - September | L | -1.9120327 | 0.4043576 | 181 | -4.7285685 | 0.0000268 |
| June - July | L | -0.6455467 | 0.4016562 | 181 | -1.6072121 | 0.3771338 |
| June - September | L | -1.6864224 | 0.3727681 | 181 | -4.5240523 | 0.0000644 |
| July - September | L | -1.0408757 | 0.3202564 | 181 | -3.2501328 | 0.0074540 |
| May - June | H | 0.0715604 | 0.5133258 | 181 | 0.1394053 | 0.9990301 |
| May - July | H | -0.1145134 | 0.4815480 | 181 | -0.2378026 | 0.9952531 |
| May - September | H | -0.9607980 | 0.4320781 | 181 | -2.2236674 | 0.1207383 |
| June - July | H | -0.1860737 | 0.4980846 | 181 | -0.3735786 | 0.9821648 |
| June - September | H | -1.0323583 | 0.4450468 | 181 | -2.3196621 | 0.0973494 |
| July - September | H | -0.8462846 | 0.4105998 | 181 | -2.0610932 | 0.1699512 |
Green Algae from transects in low and high salinity sites over the summer
Species: fucus, melanosiphon, leathesia
Notes
Mostly fucus
brown_poisson <- glmmTMB(round(browns) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
brown_NB <- update(brown_poisson, family = nbinom2)
brown_NB1 <- update(brown_poisson, family = nbinom1) #best model
brown_gaussian <- glmmTMB((browns+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
brown_tweedie <- glmmTMB(browns ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
# models with zero-inflation
brown_zipoisson <- glmmTMB(round(browns) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
brown_ziNB <- update(brown_zipoisson, family = nbinom2)
brown_ziNB1 <- update(brown_zipoisson, family = nbinom1)
Checking assumptions:
## dAIC df
## brown_NB1 0.0 11
## brown_ziNB1 1.5 12
## brown_NB 22.4 11
## brown_ziNB 24.4 12
## brown_tweedie 51.7 12
## brown_gaussian 348.9 11
## brown_zipoisson 2856.6 11
## brown_poisson 3237.8 10
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.1239, p-value = 0.8461
## alternative hypothesis: true autocorrelation is not 0
Adding a dispersion parameter for salinity_regime:
Adding dispersion parameter for site instead:
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 3.0445515 | 0.1928951 | 15.7834597 | 0.0000000 |
| monthJune | 0.1421746 | 0.2523742 | 0.5633484 | 0.5731977 |
| monthJuly | 0.6544828 | 0.2303375 | 2.8414082 | 0.0044915 |
| monthSeptember | 0.7550187 | 0.2290348 | 3.2965238 | 0.0009789 |
| salinity_regimeH | 0.2522594 | 0.3170038 | 0.7957612 | 0.4261708 |
| monthJune:salinity_regimeH | -0.4583030 | 0.4169180 | -1.0992641 | 0.2716529 |
| monthJuly:salinity_regimeH | -0.6180043 | 0.3877060 | -1.5940025 | 0.1109355 |
| monthSeptember:salinity_regimeH | -0.6116013 | 0.3798392 | -1.6101584 | 0.1073633 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 15.399269 | 3 | 0.0015054 |
| salinity_regime | 1.207230 | 1 | 0.2718814 |
| month:salinity_regime | 3.324017 | 3 | 0.3443139 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | -0.1421746 | 0.2523742 | 176 | -0.5633484 | 0.9427681 |
| May - July | L | -0.6544828 | 0.2303375 | 176 | -2.8414082 | 0.0256300 |
| May - September | L | -0.7550187 | 0.2290348 | 176 | -3.2965238 | 0.0064487 |
| June - July | L | -0.5123083 | 0.2153927 | 176 | -2.3784844 | 0.0850085 |
| June - September | L | -0.6128441 | 0.2305027 | 176 | -2.6587281 | 0.0422054 |
| July - September | L | -0.1005358 | 0.2071795 | 176 | -0.4852596 | 0.9622891 |
| May - June | H | 0.3161285 | 0.3318553 | 176 | 0.9526094 | 0.7763957 |
| May - July | H | -0.0364786 | 0.3118662 | 176 | -0.1169686 | 0.9994258 |
| May - September | H | -0.1434174 | 0.3030196 | 176 | -0.4732940 | 0.9648602 |
| June - July | H | -0.3526070 | 0.3302644 | 176 | -1.0676508 | 0.7096561 |
| June - September | H | -0.4595458 | 0.3224557 | 176 | -1.4251441 | 0.4853869 |
| July - September | H | -0.1069388 | 0.3018381 | 176 | -0.3542919 | 0.9847048 |
#### Functional Algal Groupings Grazed (sheets and filaments) vs. Non grazed (leathery, calcareous, crustose). Is there a better term for this?
Grazed Filaments, sheets Ulva, Pyropia, Polysiphonia, Diatoms, Urospora, Acrosiphonia
Non-grazed Turfs, crusts Fucus, Mastocarpus, Mastocarpus crust, Hildenbrandia, Endocladia, Melanosiphon (not sure if this is the right grouping), leathesia Endocladia and Mastocarpus both part of Lottia diet though?
grazed_algae_poisson <- glmmTMB(round(grazed) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
grazed_algae_NB <- update(grazed_algae_poisson, family = nbinom2) # Best model
grazed_algae_NB1 <- update(grazed_algae_poisson, family = nbinom1)
grazed_algae_gaussian <- glmmTMB((grazed+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
grazed_algae_tweedie <- glmmTMB(grazed ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
# models with zero-inflation
grazed_algae_zipoisson <- glmmTMB(round(grazed) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
grazed_algae_ziNB <- update(grazed_algae_zipoisson, family = nbinom2)
grazed_algae_ziNB1 <- update(grazed_algae_zipoisson, family = nbinom1)
Grazed algae: checking model assumptions
## dAIC df
## grazed_algae_NB 0.0 11
## grazed_algae_NB1 1.4 11
## grazed_algae_ziNB 2.0 12
## grazed_algae_ziNB1 2.4 12
## grazed_algae_tweedie 65.0 12
## grazed_algae_gaussian 576.2 11
## grazed_algae_zipoisson 781.4 11
## grazed_algae_poisson 1330.9 10
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.734, p-value = 0.3787
## alternative hypothesis: true autocorrelation is not 0
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 0.2750823 | 0.6570517 | 0.4186615 | 0.6754635 |
| monthJune | 0.3417898 | 0.5474157 | 0.6243698 | 0.5323847 |
| monthJuly | 1.0713042 | 0.5411664 | 1.9796208 | 0.0477462 |
| monthSeptember | 2.0280985 | 0.5326420 | 3.8076204 | 0.0001403 |
| salinity_regimeH | 1.4646335 | 0.9147464 | 1.6011361 | 0.1093468 |
| monthJune:salinity_regimeH | -0.7042222 | 0.7381904 | -0.9539844 | 0.3400915 |
| monthJuly:salinity_regimeH | -1.6986907 | 0.7425933 | -2.2875116 | 0.0221660 |
| monthSeptember:salinity_regimeH | -2.0856043 | 0.7412548 | -2.8136132 | 0.0048988 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 10.2155670 | 3 | 0.0168199 |
| salinity_regime | 0.1438865 | 1 | 0.7044475 |
| month:salinity_regime | 9.8419923 | 3 | 0.0199581 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | -0.3417898 | 0.5474157 | 181 | -0.6243698 | 0.9241470 |
| May - July | L | -1.0713042 | 0.5411664 | 181 | -1.9796208 | 0.1995108 |
| May - September | L | -2.0280985 | 0.5326420 | 181 | -3.8076204 | 0.0010930 |
| June - July | L | -0.7295144 | 0.5527433 | 181 | -1.3198070 | 0.5515370 |
| June - September | L | -1.6863087 | 0.5079056 | 181 | -3.3201222 | 0.0059442 |
| July - September | L | -0.9567942 | 0.5318073 | 181 | -1.7991371 | 0.2771120 |
| May - June | H | 0.3624323 | 0.4956890 | 181 | 0.7311688 | 0.8844714 |
| May - July | H | 0.6273865 | 0.5088229 | 181 | 1.2330153 | 0.6067162 |
| May - September | H | 0.0575059 | 0.5131515 | 181 | 0.1120641 | 0.9994949 |
| June - July | H | 0.2649541 | 0.5220343 | 181 | 0.5075416 | 0.9572101 |
| June - September | H | -0.3049265 | 0.5198372 | 181 | -0.5865807 | 0.9360315 |
| July - September | H | -0.5698806 | 0.5228355 | 181 | -1.0899806 | 0.6961113 |
non_grazed_algae_poisson <- glmmTMB(round(non_grazed) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = poisson(link = "log"))
non_grazed_algae_NB <- update(non_grazed_algae_poisson, family = nbinom2) # Best model
non_grazed_algae_NB1 <- update(non_grazed_algae_poisson, family = nbinom1)
non_grazed_algae_gaussian <- glmmTMB((non_grazed+1) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = gaussian(link = "log"))
non_grazed_algae_tweedie <- glmmTMB(non_grazed ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, family = tweedie(link = "log"))
# models with zero-inflation
non_grazed_algae_zipoisson <- glmmTMB(round(non_grazed) ~ month * salinity_regime + (1|salinity_regime/site), data = survey_data, ziformula = ~1, family = poisson(link = "log"))
non_grazed_algae_ziNB <- update(non_grazed_algae_zipoisson, family = nbinom2)
# convergence issues
non_grazed_algae_ziNB1 <- update(non_grazed_algae_zipoisson, family = nbinom1)
AICtab(non_grazed_algae_poisson, non_grazed_algae_NB, non_grazed_algae_NB1, non_grazed_algae_gaussian, non_grazed_algae_tweedie, non_grazed_algae_zipoisson, non_grazed_algae_ziNB1)
## dAIC df
## non_grazed_algae_NB1 0.0 11
## non_grazed_algae_ziNB1 1.5 12
## non_grazed_algae_NB 22.2 11
## non_grazed_algae_tweedie 48.3 12
## non_grazed_algae_gaussian 349.2 11
## non_grazed_algae_zipoisson 2852.3 11
## non_grazed_algae_poisson 3230.9 10
Checking model assumptions:
Adding a dispersion parameter for salinity regime:
Adding a dispersion formula for site instead:
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.8743, p-value = 0.8443
## alternative hypothesis: true autocorrelation is not 0
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 3.0289719 | 0.1924510 | 15.7389283 | 0.0000000 |
| monthJune | 0.1414859 | 0.2524612 | 0.5604265 | 0.5751886 |
| monthJuly | 0.6639990 | 0.2299227 | 2.8879226 | 0.0038780 |
| monthSeptember | 0.7695971 | 0.2285388 | 3.3674681 | 0.0007586 |
| salinity_regimeH | 0.2678394 | 0.3167337 | 0.8456296 | 0.3977594 |
| monthJune:salinity_regimeH | -0.4576148 | 0.4169706 | -1.0974750 | 0.2724338 |
| monthJuly:salinity_regimeH | -0.6275224 | 0.3874597 | -1.6195811 | 0.1053223 |
| monthSeptember:salinity_regimeH | -0.6261807 | 0.3795403 | -1.6498400 | 0.0989757 |
| term | statistic | df | p.value |
|---|---|---|---|
| month | 15.958681 | 3 | 0.0011563 |
| salinity_regime | 1.111575 | 1 | 0.2917398 |
| month:salinity_regime | 3.459912 | 3 | 0.3259987 |
| contrast | salinity_regime | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| May - June | L | -0.1414859 | 0.2524612 | 176 | -0.5604265 | 0.9435853 |
| May - July | L | -0.6639990 | 0.2299227 | 176 | -2.8879226 | 0.0224567 |
| May - September | L | -0.7695971 | 0.2285388 | 176 | -3.3674681 | 0.0051136 |
| June - July | L | -0.5225130 | 0.2156012 | 176 | -2.4235164 | 0.0763768 |
| June - September | L | -0.6281111 | 0.2305096 | 176 | -2.7248811 | 0.0353643 |
| July - September | L | -0.1055981 | 0.2066763 | 176 | -0.5109348 | 0.9563991 |
| May - June | H | 0.3161289 | 0.3318551 | 176 | 0.9526111 | 0.7763948 |
| May - July | H | -0.0364766 | 0.3118663 | 176 | -0.1169623 | 0.9994259 |
| May - September | H | -0.1434164 | 0.3030195 | 176 | -0.4732909 | 0.9648608 |
| June - July | H | -0.3526055 | 0.3302644 | 176 | -1.0676460 | 0.7096591 |
| June - September | H | -0.4595452 | 0.3224556 | 176 | -1.4251425 | 0.4853879 |
| July - September | H | -0.1069398 | 0.3018382 | 176 | -0.3542950 | 0.9847044 |
Abundance of non-grazed algae from transects in low and high salinity sites over the summer
## Field Experiment Data Notes
Only using the control and exclusion plots.
Not modelling limpets b/c they will be different due to treatment.
Do we just assume the herbivores that were removed from the exclusion plots were not having an impact on the community? i.e. they weren’t there for long enough?
Waiting to hear back from Theraesa about:
Additional Limpets in ring (Removed) two of the control plots show 34 and 1 paradigitalis in this column but 0s under the paradigitalis column.
Total Limpets (Non-Pelta removed from rings) some of the control plots have other limpet species as removed? Littorina Removed from Ring (#) some of the control plots have littorines removed?
Species: Balanus, Chthamalus, limpets, littorina, mytilus, amphipods diatoms, pyropia, fucus, ulva,
Check and see whether the exclusion plots had less limpets than the controls. We expect the interaction to be significant as well. I’m using data from May, June and July for this as previous months will have had an impact on July’s communities as well.
## dAIC df
## exp_limpets_NB 0.0 7
## exp_limpets_NB1 10.4 7
## exp_limpets_tweedie 29.4 8
## exp_limpets_zipoisson 81.0 7
## exp_limpets_poisson 391.0 6
## exp_limpets_gaussian 1079.8 7
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -3.2938817 | 0.9720340 | -3.3886488 | 0.0007024 |
| treatmentE | 0.6051422 | 0.9769437 | 0.6194238 | 0.5356372 |
| regionHigh | 3.2174139 | 1.2048968 | 2.6702817 | 0.0075788 |
| treatmentE:regionHigh | -1.8579827 | 1.1802315 | -1.5742529 | 0.1154290 |
| term | statistic | df | p.value |
|---|---|---|---|
| treatment | 1.519823 | 1 | 0.2176463 |
| region | 4.699369 | 1 | 0.0301737 |
| treatment:region | 2.478272 | 1 | 0.1154290 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| C - E | Low | -0.6051422 | 0.9769437 | 245 | -0.6194238 | 0.5362128 |
| C - E | High | 1.2528406 | 0.6580842 | 245 | 1.9037694 | 0.0581120 |
Question How were barnacles counted? How can we have non-integers?
exp_balanus_poisson <- glmmTMB(round(balanus_no) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_balanus_NB <- update(exp_balanus_poisson, family = nbinom2)
exp_balanus_NB1 <- update(exp_balanus_poisson, family = nbinom1)
exp_balanus_gaussian <- glmmTMB((balanus_no+1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_balanus_tweedie <- glmmTMB(balanus_no ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log")) # Best model
# models with zero-inflation
exp_balanus_zipoisson <- glmmTMB(round(balanus_no) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_balanus_ziNB <- update(exp_balanus_zipoisson, family = nbinom2) # convergence issues
exp_balanus_ziNB1 <- update(exp_balanus_zipoisson, family = nbinom1) # convergence issues
Checking model assumptions:
## dAIC df
## exp_balanus_tweedie 0.0 8
## exp_balanus_NB1 2.4 7
## exp_balanus_gaussian 10.0 7
## exp_balanus_NB 20.4 7
## exp_balanus_poisson 13064.6 6
## exp_balanus_zipoisson 13066.6 7
Adding a dispersion formula for trmt*region:
Trying to fit dispersion formula to treatment*site but it’s causing convergence issues.
What if I don’t include site as a random effect, but as a fixed effect instead?
Adding a dispersion parameter for treatment*site
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 5.1572444 | 0.3390990 | 15.2086678 | 0.0000000 |
| treatmentexclusion | 0.6170414 | 0.3796994 | 1.6250790 | 0.1041457 |
| siteRope Site 2 | 2.0630250 | 0.3410859 | 6.0484033 | 0.0000000 |
| siteHorseshoe Bay | 0.8385693 | 0.4697243 | 1.7852373 | 0.0742228 |
| siteRuckle | 1.7846913 | 0.4240018 | 4.2091596 | 0.0000256 |
| siteEagle Cove | 2.3308623 | 0.3667312 | 6.3557795 | 0.0000000 |
| siteHailstorm 2 | 1.2292992 | 0.4255567 | 2.8886850 | 0.0038686 |
| treatmentexclusion:siteRope Site 2 | -0.6117046 | 0.3827412 | -1.5982199 | 0.1099941 |
| treatmentexclusion:siteHorseshoe Bay | -0.0399731 | 0.5254315 | -0.0760767 | 0.9393581 |
| treatmentexclusion:siteRuckle | -0.8247594 | 0.4626999 | -1.7824933 | 0.0746688 |
| treatmentexclusion:siteEagle Cove | -0.8251589 | 0.4296878 | -1.9203683 | 0.0548114 |
| treatmentexclusion:siteHailstorm 2 | -0.2195158 | 0.4683627 | -0.4686877 | 0.6392928 |
| term | statistic | df | p.value |
|---|---|---|---|
| (Intercept) | 231.303576 | 1 | 0.0000000 |
| treatment | 2.640882 | 1 | 0.1041457 |
| site | 65.318524 | 5 | 0.0000000 |
| treatment:site | 8.834319 | 5 | 0.1158579 |
| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| Lions Bay 2 - Rope Site 2 | -1.7571727 | 0.1913706 | 59 | -9.182041 | 0.0000000 |
| Lions Bay 2 - Horseshoe Bay | -0.8185828 | 0.2627158 | 59 | -3.115850 | 0.0321635 |
| Lions Bay 2 - Ruckle | -1.3723116 | 0.2313499 | 59 | -5.931757 | 0.0000025 |
| Lions Bay 2 - Eagle Cove | -1.9182829 | 0.2148439 | 59 | -8.928728 | 0.0000000 |
| Lions Bay 2 - Hailstorm 2 | -1.1195413 | 0.2341813 | 59 | -4.780660 | 0.0001691 |
| Rope Site 2 - Horseshoe Bay | 0.9385899 | 0.1831843 | 59 | 5.123747 | 0.0000495 |
| Rope Site 2 - Ruckle | 0.3848610 | 0.1343864 | 59 | 2.863840 | 0.0612430 |
| Rope Site 2 - Eagle Cove | -0.1611102 | 0.1034158 | 59 | -1.557888 | 0.6288773 |
| Rope Site 2 - Hailstorm 2 | 0.6376314 | 0.1392042 | 59 | 4.580547 | 0.0003403 |
| Horseshoe Bay - Ruckle | -0.5537288 | 0.2246254 | 59 | -2.465121 | 0.1514418 |
| Horseshoe Bay - Eagle Cove | -1.0997001 | 0.2075854 | 59 | -5.297580 | 0.0000263 |
| Horseshoe Bay - Hailstorm 2 | -0.3009585 | 0.2275405 | 59 | -1.322659 | 0.7713864 |
| Ruckle - Eagle Cove | -0.5459712 | 0.1661171 | 59 | -3.286664 | 0.0202130 |
| Ruckle - Hailstorm 2 | 0.2527703 | 0.1904675 | 59 | 1.327104 | 0.7689080 |
| Eagle Cove - Hailstorm 2 | 0.7987416 | 0.1700382 | 59 | 4.697423 | 0.0002266 |
Balanus (percent cover) in control vs. exclusion plots in low and high salinity sites
exp_chthamalus_poisson <- glmmTMB(round(chthamalus_no) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_chthamalus_NB <- update(exp_chthamalus_poisson, family = nbinom2)
exp_chthamalus_NB1 <- update(exp_chthamalus_poisson, family = nbinom1)
exp_chthamalus_gaussian <- glmmTMB((chthamalus_no+1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_chthamalus_tweedie <- glmmTMB(chthamalus_no ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log"))
# models with zero-inflation
exp_chthamalus_zipoisson <- glmmTMB(round(chthamalus_no) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_chthamalus_ziNB <- update(exp_chthamalus_zipoisson, family = nbinom2)
exp_chthamalus_ziNB1 <- update(exp_chthamalus_zipoisson, family = nbinom1) # best model
Model Selection and checking assumptions:
## dAIC df
## exp_chthamalus_ziNB1 0.0 8
## exp_chthamalus_tweedie 3.2 8
## exp_chthamalus_NB1 5.7 7
## exp_chthamalus_ziNB 19.0 8
## exp_chthamalus_NB 27.8 7
## exp_chthamalus_gaussian 402.8 7
## exp_chthamalus_zipoisson 2222.5 7
## exp_chthamalus_poisson 3612.3 6
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 2.3549717 | 1.0679766 | 2.2050780 | 0.0274486 |
| treatmentexclusion | 0.0791074 | 0.3619432 | 0.2185629 | 0.8269905 |
| regionHigh | 2.3623231 | 1.4142626 | 1.6703567 | 0.0948488 |
| treatmentexclusion:regionHigh | -2.2124492 | 0.5343787 | -4.1402269 | 0.0000347 |
| term | statistic | df | p.value |
|---|---|---|---|
| treatment | 12.613247 | 1 | 0.0003830 |
| region | 1.441596 | 1 | 0.2298812 |
| treatment:region | 17.141479 | 1 | 0.0000347 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | -0.0791074 | 0.3619432 | 76 | -0.2185629 | 0.8275764 |
| control - exclusion | High | 2.1333418 | 0.3913361 | 76 | 5.4514317 | 0.0000006 |
Chthamalus (percent cover) in control vs. exclusion plots in low and high salinity sites
## # A tibble: 6 × 2
## site n
## <fct> <int>
## 1 Lions Bay 2 14
## 2 Rope Site 2 14
## 3 Horseshoe Bay 14
## 4 Ruckle 11
## 5 Eagle Cove 8
## 6 Hailstorm 2 12
Fit models to all algae pooled together:
exp_algae_poisson <- glmmTMB(round(algae) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_algae_NB <- update(exp_algae_poisson, family = nbinom2)
exp_algae_NB1 <- update(exp_algae_poisson, family = nbinom1)
exp_algae_gaussian <- glmmTMB((algae + 1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_algae_tweedie <- glmmTMB(algae ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log"))
# models with zero-inflation
exp_algae_zipoisson <- glmmTMB(round(algae) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_algae_ziNB <- update(exp_algae_zipoisson, family = nbinom2)
exp_algae_ziNB1 <- update(exp_algae_zipoisson, family = nbinom1) #best model
Checking model assumptions:
## dAIC df
## exp_algae_ziNB1 0.0 8
## exp_algae_NB1 5.8 7
## exp_algae_tweedie 15.8 8
## exp_algae_ziNB 34.7 8
## exp_algae_NB 41.9 7
## exp_algae_gaussian 76.5 7
## exp_algae_zipoisson 764.4 7
## exp_algae_poisson 1117.8 6
Adding a dispersion parameter for region:
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 4.045601 | 0.2073427 | 19.5116589 | 0.0000000 |
| treatmentexclusion | -0.115691 | 0.1551350 | -0.7457439 | 0.4558221 |
| regionHigh | -2.179333 | 0.4570558 | -4.7681986 | 0.0000019 |
| treatmentexclusion:regionHigh | 2.274326 | 0.4043492 | 5.6246579 | 0.0000000 |
| term | statistic | df | p.value |
|---|---|---|---|
| treatment | 2.3728287 | 1 | 0.1234632 |
| region | 0.8474256 | 1 | 0.3572816 |
| treatment:region | 31.6367762 | 1 | 0.0000000 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | 0.115691 | 0.1551350 | 75 | 0.7457439 | 0.4581531 |
| control - exclusion | High | -2.158635 | 0.3731297 | 75 | -5.7852129 | 0.0000002 |
Algae in control vs. exclusion plots in low and high salinity sites
Almost all Mastocarpus crust with a very small amount of Pyropia
exp_red_poisson <- glmmTMB(round(reds) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_red_NB <- update(exp_red_poisson, family = nbinom2)
exp_red_NB1 <- update(exp_red_poisson, family = nbinom1)
exp_red_gaussian <- glmmTMB((reds + 1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_red_tweedie <- glmmTMB(reds ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log")) # best model
# models with zero-inflation
exp_red_zipoisson <- glmmTMB(round(reds) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_red_ziNB <- update(exp_red_zipoisson, family = nbinom2)
exp_red_ziNB1 <- update(exp_red_zipoisson, family = nbinom1)
## dAIC df
## exp_red_tweedie 0.0 8
## exp_red_zipoisson 4.4 7
## exp_red_ziNB1 4.9 8
## exp_red_ziNB 6.1 8
## exp_red_NB1 9.8 7
## exp_red_NB 16.7 7
## exp_red_poisson 73.3 6
## exp_red_gaussian 172.5 7
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.0234567 | 0.7880334 | -0.0297661 | 0.9762536 |
| treatmentexclusion | -0.3085752 | 0.4686554 | -0.6584266 | 0.5102640 |
| regionHigh | -0.4996745 | 1.0881564 | -0.4591937 | 0.6460951 |
| treatmentexclusion:regionHigh | -0.2724478 | 0.8529393 | -0.3194223 | 0.7494063 |
| term | statistic | df | p.value |
|---|---|---|---|
| treatment | 0.9961414 | 1 | 0.3182460 |
| region | 0.3475700 | 1 | 0.5554919 |
| treatment:region | 0.1020306 | 1 | 0.7494063 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | 0.3085752 | 0.4686554 | 76 | 0.6584266 | 0.5122538 |
| control - exclusion | High | 0.5810230 | 0.7126650 | 76 | 0.8152820 | 0.4174586 |
Red Algae in control vs. limpet exclusion plots in high and low salinity sites
Only green algae species is Ulva
Only region is significant, because the confidence intervals are very large.
How many plots, and in which sites, have green algae?
## # A tibble: 6 × 2
## site n
## <fct> <int>
## 1 Lions Bay 2 14
## 2 Rope Site 2 14
## 3 Horseshoe Bay 12
## 4 Ruckle 8
## 5 Eagle Cove 6
## 6 Hailstorm 2 5
exp_green_poisson <- glmmTMB(round(greens) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_green_NB <- update(exp_green_poisson, family = nbinom2)
exp_green_NB1 <- update(exp_green_poisson, family = nbinom1)
exp_green_gaussian <- glmmTMB((greens+1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_green_tweedie <- glmmTMB(greens ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log"))
# models with zero-inflation
exp_green_zipoisson <- glmmTMB(round(greens) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_green_ziNB <- update(exp_green_zipoisson, family = nbinom2)
exp_green_ziNB1 <- update(exp_green_zipoisson, family = nbinom1) # best model
## dAIC df
## exp_green_ziNB1 0.0 8
## exp_green_tweedie 11.4 8
## exp_green_NB1 14.1 7
## exp_green_ziNB 43.3 8
## exp_green_NB 79.1 7
## exp_green_gaussian 166.8 7
## exp_green_zipoisson 539.5 7
## exp_green_poisson 1310.0 6
Adding a dispersion formula for region:
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 4.0687806 | 0.1415476 | 28.7449540 | 0.0000000 |
| treatmentexclusion | -0.1154469 | 0.1417591 | -0.8143877 | 0.4154229 |
| regionHigh | -3.5833048 | 0.8293180 | -4.3207852 | 0.0000155 |
| treatmentexclusion:regionHigh | 3.5105808 | 0.8378930 | 4.1897723 | 0.0000279 |
| term | statistic | df | p.value |
|---|---|---|---|
| (Intercept) | 826.2723828 | 1 | 0.0000000 |
| treatment | 0.6632273 | 1 | 0.4154229 |
| region | 18.6691845 | 1 | 0.0000155 |
| treatment:region | 17.5541920 | 1 | 0.0000279 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | 0.1154469 | 0.1417591 | 75 | 0.8143877 | 0.4180012 |
| control - exclusion | High | -3.3951339 | 0.8248749 | 75 | -4.1159380 | 0.0000980 |
Only brown algae species is Fucus
How many plots, and in which sites, have brown algae?
No brown algae at Ruckle Park
## # A tibble: 5 × 2
## site n
## <fct> <int>
## 1 Lions Bay 2 3
## 2 Rope Site 2 2
## 3 Horseshoe Bay 3
## 4 Eagle Cove 2
## 5 Hailstorm 2 1
exp_brown_poisson <- glmmTMB(round(browns) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_brown_NB <- update(exp_brown_poisson, family = nbinom2)
exp_brown_NB1 <- update(exp_brown_poisson, family = nbinom1)
exp_brown_gaussian <- glmmTMB((browns + 1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_brown_tweedie <- glmmTMB(browns ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log")) # best model
# models with zero-inflation
exp_brown_zipoisson <- glmmTMB(round(browns) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_brown_ziNB <- update(exp_brown_zipoisson, family = nbinom2) #convergence issues
exp_brown_ziNB1 <- update(exp_brown_zipoisson, family = nbinom1)
Checking assumptions:
## dAIC df
## exp_brown_tweedie 0.0 8
## exp_brown_zipoisson 7.3 7
## exp_brown_ziNB1 9.3 8
## exp_brown_NB1 10.8 7
## exp_brown_NB 11.0 7
## exp_brown_poisson 13.5 6
## exp_brown_gaussian 25.9 7
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.8934872 | 4.104650e-01 | -2.1767684 | 0.0294979 |
| treatmentexclusion | -1.0116009 | 7.677946e-01 | -1.3175411 | 0.1876573 |
| regionHigh | -1.2992830 | 8.454432e-01 | -1.5368069 | 0.1243406 |
| treatmentexclusion:regionHigh | -21.3472632 | 1.932017e+04 | -0.0011049 | 0.9991184 |
| term | statistic | df | p.value |
|---|---|---|---|
| (Intercept) | 4.7383208 | 1 | 0.0294979 |
| treatment | 1.7359146 | 1 | 0.1876573 |
| region | 2.3617753 | 1 | 0.1243406 |
| treatment:region | 0.0000012 | 1 | 0.9991184 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | 1.011601 | 7.677946e-01 | 76 | 1.3175411 | 0.1916155 |
| control - exclusion | High | 22.358864 | 1.932017e+04 | 76 | 0.0011573 | 0.9990797 |
Brown Algae in control vs. limpet exclusion plots in high and low salinity sites
Only 2 sites have diatoms so not doing anything with this at the moment.
## # A tibble: 2 × 2
## site n
## <fct> <int>
## 1 Ruckle 6
## 2 Hailstorm 2 2
exp_grazed_poisson <- glmmTMB(round(grazed) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_grazed_NB <- update(exp_grazed_poisson, family = nbinom2)
exp_grazed_NB1 <- update(exp_grazed_poisson, family = nbinom1)
exp_grazed_gaussian <- glmmTMB((grazed+1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_grazed_tweedie <- glmmTMB(grazed ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log"))
# models with zero-inflation
exp_grazed_zipoisson <- glmmTMB(round(grazed) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_grazed_ziNB <- update(exp_grazed_zipoisson, family = nbinom2)
exp_grazed_ziNB1 <- update(exp_grazed_zipoisson, family = nbinom1) # best model
Checking assumptions:
## dAIC df
## exp_grazed_ziNB1 0.0 8
## exp_grazed_tweedie 11.0 8
## exp_grazed_NB1 18.7 7
## exp_grazed_ziNB 41.7 8
## exp_grazed_NB 81.0 7
## exp_grazed_gaussian 135.6 7
## exp_grazed_zipoisson 542.0 7
## exp_grazed_poisson 1291.1 6
Adding a dispersion formula for treatment*site
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 4.0771477 | 0.1244534 | 32.7604248 | 0.0000000 |
| treatmentexclusion | -0.1053661 | 0.1454185 | -0.7245718 | 0.4687148 |
| regionHigh | -1.2498128 | 1.0176401 | -1.2281482 | 0.2193914 |
| treatmentexclusion:regionHigh | 1.6313215 | 1.0190669 | 1.6007992 | 0.1094214 |
| term | statistic | df | p.value |
|---|---|---|---|
| (Intercept) | 1073.2454346 | 1 | 0.0000000 |
| treatment | 0.5250043 | 1 | 0.4687148 |
| region | 1.5083479 | 1 | 0.2193914 |
| treatment:region | 2.5625582 | 1 | 0.1094214 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | 0.1053661 | 0.1454185 | 65 | 0.7245718 | 0.4713150 |
| control - exclusion | High | -1.5259554 | 1.0092699 | 65 | -1.5119399 | 0.1353961 |
exp_nongrazed_poisson <- glmmTMB(round(nongrazed) ~ treatment * region + (1|region/site), data = exp_data, family = poisson(link = "log"))
exp_nongrazed_NB <- update(exp_nongrazed_poisson, family = nbinom2)
exp_nongrazed_NB1 <- update(exp_nongrazed_poisson, family = nbinom1)
exp_nongrazed_gaussian <- glmmTMB((nongrazed+1) ~ treatment * region + (1|region/site), data = exp_data, family = gaussian(link = "log"))
exp_nongrazed_tweedie <- glmmTMB(nongrazed ~ treatment * region + (1|region/site), data = exp_data, family = tweedie(link = "log")) # best model
# models with zero-inflation
exp_nongrazed_zipoisson <- glmmTMB(round(nongrazed) ~ treatment * region + (1|region/site), data = exp_data, ziformula = ~1, family = poisson(link = "log"))
exp_nongrazed_ziNB <- update(exp_nongrazed_zipoisson, family = nbinom2)
exp_nongrazed_ziNB1 <- update(exp_nongrazed_zipoisson, family = nbinom1)
Checking model assumptions:
## dAIC df
## exp_nongrazed_tweedie 0.0 8
## exp_nongrazed_ziNB1 6.8 8
## exp_nongrazed_NB1 7.2 7
## exp_nongrazed_ziNB 13.4 8
## exp_nongrazed_NB 14.9 7
## exp_nongrazed_zipoisson 15.7 7
## exp_nongrazed_poisson 54.4 6
## exp_nongrazed_gaussian 142.6 7
| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 0.4324789 | 0.5585112 | 0.7743425 | 0.4387283 |
| treatmentexclusion | -0.4258574 | 0.3534379 | -1.2049000 | 0.2282419 |
| regionHigh | -0.7537305 | 0.8135738 | -0.9264440 | 0.3542153 |
| treatmentexclusion:regionHigh | -0.4247930 | 0.6973658 | -0.6091394 | 0.5424320 |
| term | statistic | df | p.value |
|---|---|---|---|
| (Intercept) | 0.5996063 | 1 | 0.4387283 |
| treatment | 1.4517840 | 1 | 0.2282419 |
| region | 0.8582985 | 1 | 0.3542153 |
| treatment:region | 0.3710508 | 1 | 0.5424320 |
| contrast | region | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| control - exclusion | Low | 0.4258574 | 0.3534379 | 76 | 1.204900 | 0.2319805 |
| control - exclusion | High | 0.8506504 | 0.6012904 | 76 | 1.414708 | 0.1612383 |
| species | distribution_family | zero_inflated | dispersion_formula | month | region | month_region |
|---|---|---|---|---|---|---|
| limpets | negative binomial 1 | yes | n/a | yes | yes | yes |
| littorines | negative binomial 1 | no | salinity*month | no | no | no |
| balanus | negative binomial 1 | no | salinity*month | yes | yes | yes |
| chthamalus | negative binomial 2 | no | salinity | yes | yes | yes |
| red algae | negative binomial 1 | no | site | no | yes | no |
| green algae | negative binomial 1 | no | n/a | yes | yes | yes |
| brown algae | negative binomial 1 | no | site | yes | no | no |
| grazed | negative binomial 2 | no | n/a | yes | no | yes |
| nongrazed | negative binomial 1 | no | site | yes | no | no |
| species | distribution_family | zero_inflated | dispersion_formula | treatment | region | treatment_region |
|---|---|---|---|---|---|---|
| balanus | tweedie | no | treatment*site | yes | yes (site) | yes |
| chthamalus | negative binomial 1 | yes | n/a | yes | yes | yes |
| all algae | negative binomial 1 | yes | region | yes | yes | yes |
| red algae | tweedie | no | n/a | no | no | no |
| green algae | negative binomial 1 | yes | region | no | yes | no |
| brown algae | tweedie | no | n/a | no | no | no |
| grazed | negative binomial 1 | yes | treatment*site | no | yes | no |
| nongrazed | tweedie | no | n/a | no | no | no |